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ABSTRACT 

While recent observational progress is converging on the detection of compact regions of thermal emis¬ 
sion due to embedded protoplanets, further theoretical predictions are needed to understand the response 
of a protoplanetary disk to the planet formation radiative feedback. This is particularly important to 
make predictions for the observability of circumplanetary regions. In this work we use 2D hydrodynamical 
simulations to examine the evolution of a viscous protoplanetary disk in which a luminous Jupiter-mass 
planet is embedded. We use an energy equation which includes the radiative heating of the planet as an 
additional mechanism for planet formation feedback. Several models are computed for planet luminosities 
ranging from 10 -5 to 10 -3 Solar luminosities. We find that the planet radiative feedback enhances the 
disk’s accretion rate at the planet’s orbital radius, producing a hotter and more luminous environement 
around the planet, independently of the prescription used to model the disk’s turbulent viscosity. We also 
estimate the thermal signature of the planet feedback for our range of planet luminosities, finding that the 
emitted spectrum of a purely active disk, without passive heating, is appreciably modified in the infrared. 

We simulate the protoplanetary disk around HD 100546 where a planet companion is located at about 
68 au from the star. Assuming the planet mass is 5 Jupiter masses and its luminosity is ^ 2.5 x 10 -4 L©, 
we find that the radiative feedback of the planet increases the luminosity of its ^ 5 au circumplanetary 
disk from 10 -5 L 0 (without feedback) to 10 _3 L o , corresponding to an emission of ^ 1 mJy in L' band 
after radiative transfer calculations, a value that is in good agreement with HD 100546b observations. 

Subject headings: planet-disk interactions - protoplanetary disks - accretion, accretion disks - planetary 
systems - hydrodynamics - methods: numerical 


1. INTRODUCTION 

Recent observational progress allows the detailed study 
of planet formation feedback. The ALMA facility has 
opened the resolved study of accretion kinematics in 
pr otop la netary gaps. Fo r instance, Casa ssus et al. 
(|Casassus et al.l (|2012l b ICasassns et all (|2013l )) find that 
the dust gap in the disk around the star HD 142527 
shows a disrupted outer disk suggestive of on-going dy¬ 
namical clearing, and contains residual gas whose kine¬ 
matics are consistent with accretion across the dust gap. 
Dramatic advances in high-contrast imaging techniques 
have allowed the li kely de tection of embedded protoplan¬ 
ets. IQuanz et al.l (|2013af ) found a compact but resolved 
3.8/im ( L') source a t ~ 68 an from HD 100546 (indepen¬ 
dently confirmed by I Currie et ah (|2014f ), which could be 
interpreted as on-going accretion o nto a compact body . 
In another example of resolved data, IQuanz et al. 62013b ) 
present polarized light images of HD 169142 resolving 
features in its protoplanetary disk that could be inter¬ 
preted as a gap ind u ced b y for ming protop l anets . In¬ 
deed, IRnggiani et al.l (|2014f ) and IBiller et~ahl (|2014f ) find 
an L' point source, within this gap, at a separation of 
~ 22 au. There are indications that this L' compact signal 
is not photospheric, and that it is somehow connected to 
the protoplanetary accretion luminosity. Similar findings 
have been reported for the com pact Hg sig nal found at 
12 au from HD 142527 by I Close et al.l ( 2014 b which coin¬ 
cides with a relatively bright L' signal. That source would 
reach the stellar mass regime for photospheric emission 
(Bille r et al.l 1201 4b but has been re cently resolved to b e 
extended and polarized in Y band (Rodig as et al.l [2014b 
suggesting that this companion to HD 142527 is proba¬ 


bly a substellar object with a remarkably strong thermal 
luminosity in Z/, that is somehow connected to accretion. 

Hydrodynamical simulations can be used to inform the 
interpretation of the data on planet forming systems, by 
calculating the radiative emissions of a disk with an em¬ 
bedded planet t hat has created a gap. T his was the ap¬ 
proach used by I Wolf fe D’Angelo! (2005), who modelled 
in two dimensions the gravitational response of a disk 
to the presence of an embedded planet, under the as¬ 
sumption of a constant aspect ratio for the disk. As a 
post- processing step, once the simulation reached steady 
state, I Wolf fe D’Angelol ( 2005 ) assumed a luminosity for 
the planet and calculated the radiative response of the 
disk. They concluded that a hot circumplanetary region 
could eventually be detected by ALMA. 

In this paper we also use 2D hydrodynamical simula¬ 
tions to follow the dynamics of a disk with an embedded 
planet. However, instead of assuming a temperature pro¬ 
file for the disk, we use a non-stationary energy equation 
that includes the radiative feedback of planet formation 
and a temperature-dependent black body cooling for the 
disk. 

The structure of the paper is as follows: In $2] we present 
the model including the main assumptions, the physical 
conditions of the disk, as well as the numerical setup and a 
description of the code. Our results and main conclusions 
of the evolution of the density profile, temperature, and 
the spectral signature of the disk are presented in 0 with 
a short discussion in 0 We summarize our findings in £0 


2. THE MODEL 
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We are interested in the evolution of a gaseous proto¬ 
planetary disk in which a luminous Jupiter-mass planet is 
embedded. In our simulations we use an energy equation 
that includes radiative cooling, and both viscous heating 
and heating due to the planet’s luminosity. We follow the 
evolution of the protoplanet ary disk for about 10 4 years, 
assuming the planet is already formed at the beginning 
of the simulations. This is a short period of time com¬ 
pared wit h both the lifetim e of pr otoplanetary disks 
Myr e.g., Williams & Cieza (2011)), and t he time-scales 
over which planet lumin osities should vary (M arlev et al.l 
(|2007[ ) . iMordasinil (j2013h ) . thus justifying the use of a con¬ 
stant planet luminosity in the simulations. 

Our simulations use the public two- dimensional hydro¬ 
dynamics code FARGO- AE0 (IBaruteau fc Masse! 112008 ) 
which is dedicated to planet-disk interactions. It is a 
staggered mesh code that solves the Navier-Stokes, con¬ 
tinuity and energy equations on a polar grid. It is based 
on an Eulerian formalism using a finite difference method 
of second order, according to the Ivan LeeJ (|1977[ ) upwind 
algorit hm. Details of the code can b e found in [Masset 
( 2000 ) and IBaruteau fc Massetl (2008). In FARGO-AD’s 
public version, the energy equation includes viscous heat¬ 
ing and a simple temperature relaxation to reach thermo¬ 
dynamical equilibrium over some (user-defined) character¬ 
istic timescale. 

The present work features two main changes to the en¬ 
ergy equation. One is the inclusion of a radiative cooling 
function, based on the assumption that the disk radiates 
locally as a blackbody. The second, and most important 
feature, is the implementation of a heating source term as¬ 
sociated to the planet. We assume that the protoplanet 
has an intrinsic constant luminosity, therefore, it injects 
energy into the disk at a constant rate. We only take into 
account the gravitational potentials of the star and of the 
planet, the disk’s self-gravity is neglected. 

2.1. Code units and Initial Setup 

We set the mass of the central star (M*) and the planet’s 
orbital radius (r p ) as the code’s units of mass and length, 
respectively. The code’s unit of time (to) is the planet’s 
orbital period divided by 27r, that is to = (GM*/r 3 ) -1 / 2 . 
The gravitational constant G = 1 in code units. The code’s 
unit of temperature is GM*tmi p /(kBr p ), with fi the mean 
molecular weight of the gas (/x = 2.35 in all our simu¬ 
lations), m p the proton mass and ks the Boltzmann con¬ 
stant. Unless otherwise noted (see §EU), we adopt a solar- 
mass star (M*. = M 0 ) and a planet at r p = 10 au. 

We use cylindrical coordinates (r, 0). The computational 
domain extends from r = 1 to 50 au over n r = 400 equally 
spaced radial rings. It covers the full 27 t extent in azimuth 
over 77,0 = 800 equally spaced azimuthal sectors. Tests 
with higher grid resolutions,n r x = 512 x 1536, were 
performed to check the convergence of our results. We use 
an open inner and outer boundary conditions, meaning 
that the material is allowed to outflow at the disk edges. 

The initial density profile scales with r -1 : 

S(r) = £ A (1) 

r 

where Eo = 2.56 x lO _5 M 0 /au 2 . This initial disk mass is 
thus 10 -3 M q . The disk’s aspect ratio h = c s /v k, with c s 
the isothermal sound speed and vk the Keplerian velocity, 

1 http://fargo.in2p3.fr/spip.php?rubrique9 


initially equals 0.05, uniformly. The disk’s initial temper¬ 
ature therefore decreases in r -1 and ~ 63 K at 10 au for 
our fiducial primary mass (M* = Mq). 

We fix the planet-to-primary mass ratio (q) to q = 10 -3 , 
so the planet has a Jovian mass for a Solar-mass star. 
The planet is held on a fixed circular orbit (it does not 
migrate through the disk). To avoid a violent response of 
the disk to the planet’s gravitational potential initially, the 
planet mass is increased gradually over the first five orbits 
according to 

M(t) = M p sin 2 (7rf/10T p ), (2) 

where T p is the planet’s orbital period. We adopt three 
values for th e planet luminosity: 10 —5 , 10 -4 , and 10 -3 L 0 
(see § I2.3.2p . For the viscosity prescription, in our fidu¬ 
cial model we use an alpha disk model lShakura fc Sunvaevl 
(|l973l ) setting a = 4x 10 -3 (but see § 13. 4p . 

2.2. The energy equation 

The energy equation satisfied by the thermal energy den¬ 
sity e reads (e.g., citeD’Angelo-et-al-2003) 

^ + ^-(eV) = -W-V + Q+-Q-, (3) 

where it is the gas velocity, P the pressure, Q + the heat¬ 
ing rate per unit area, and Q~ the radiative cooling rate 
per unit area. To close the system of equations, an ideal 
equation of state is used, 

P = ETR, (4) 

with T the gas temperature and R = ks/iam p . The ther¬ 
mal energy density is related to the temperature through 

e=Er (A)' (5) 

where 7 denotes the adiabatic index, which we fix to 7 = 
1.4 (a typical value for a diatomic gas). Eq. (j3j) can be 
recast as 

^7 + ^ • (eV) = -(7 - l)e\^ • V + 
ot 

Qt + Qt~Q~, ( 6 ) 

where Qf is the viscous heating rate, Q p the flux of radia¬ 
tive energy received from the planet (feedback), and Q~ 
corresponds to the radiative cooling rate of the disk. These 
source terms are detailed in the next section. 

In a more realistic situation, a thermal diffusion flux 
term should be added on the right side of Eq. (J6]) . But, 
since the outer temperature of the CPD matches that of 
the protoplanet ary disk, we expect no strong heat diffusion 
to occur. 

2.3. Sources of heating and cooling 
2.3.1. Viscous dissipation 

The viscous heating rate implemented in FARGO-AD, 
Qv , has two contributions. The first one arises from the 
shear kinematic viscosity z/, 

Qiear = [ T r,r + 2r r,0 + T IA + “"g - ( 7 ) 
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where the r a $ are the components of the viscous stress 
tensor: 


t r , r — 

= 2z^Xl 


dv, 

d 

d 




0_ (V£\ 1 dvr 
dr V r / r dd> 


lto± + ±-li} .T) 

r o<fi r 3 


(8) 


Note that for a Keplerian disk (i.e., v r = 0, = vQk with 

VLk the Keplerian angular frequency), r r?r = = 0 and 

the shear viscous heating rate reduces to Q+ = = 

| vYM 2 k . 

The second contribution to the viscous heating rate is 
through the use of a von Neu mann-Richtmyer artificia l 
bulk viscosity, as described in IStone fe Normanl ([1992), 
where the coefficient C 2 is taken equal to 1.4 (C 2 mea¬ 
sures the number of zones over which a shock is spread 
over by the artificial viscosity). 


2.3.2. Planet feedback 

The key feature in this work is the inclusion of the planet 
feedback (i.e., the energy flux received from the proto¬ 
planet by the gas disk) during its evolution. In this initial 
study, we use a simplified model for the planet luminosity, 
which is powered by its mass build up during its forma¬ 
tion and evolution. Assuming free-fall, the typical power 
dissipated during the formation of a Jupiter-like planet is 


Ej 


GMj 

Rjtp ’ 


(9) 


where Mj and Rj are the mass and radius of Jupiter, re¬ 
spectively, and t p is a characteristic timescale for its for¬ 
mation. 

A characteristic timescale obtained from ground-based 
and Spitzer -based infrared (IR) surveys of young stellar 
clusters, which trace the evolution of primordial protoplan¬ 
etary disks, s uggest a formation period of about ^ 3 x 10 6 
yr (IMamajekl 2009 b Then, the emitted energy (Eq. [9j) 
reads Ej ~ 1 x 10 -4 L©, which also agrees with h ot ac¬ 
cretion shock structure formation models (e.g.. IMordasinil 
( 2013 )). To account for the large uncertainties in the accre¬ 
tion process of planetesimals, we adopt planet luminosities 
Lp ranging from 10 _1 Fj to lOEj (i.e., 10 -5 - 10 -3 L 0 ). 

In our model, the planet has already reached its final 
mass at the beginning of the simulation, and no longer 
grows. In that sense the source of the planet luminos¬ 
ity should be thought of as a “post-formation” slow con¬ 
traction of the planet, rather than accretion luminosity. 
However, it is still possible t hat with opacities larger than 
those commonly used (e.g., IBell fc Lml (| 19941 )) the accre¬ 
tion shock energy takes longer to be released, and therefore 
the total planet luminosity might include contributions of 
both aforementioned sources. Regarding the luminosities, 
it should be noted that the high values we are using (i.e., 
L p = 10 -3 L©) are expected for large opacities under the 
ass umption of a hot accretion formation model (jMordasinil 
120131 ). This is different with cold accretion initial condi¬ 
tions, which would produce smaller post-formation planet 
luminosities (< 10 -4 L 0 ). As current models of giant 
planet formation cannot distinguish between them, we as¬ 
sume hot accretion and study the effect of high planet 


luminosities. 

In all our simulations, the planet’s circumplanetary disk 
(CPD) is optically thick. The thermal energy released by 
the planet is therefore chosen to be distributed isotrop¬ 
ically within the planet’s CPD, the size of_which is de¬ 
noted by Rcpd- Following ICrida et all ( 2009 ). we adopt 

Rcpd = 0.6-Rhui, where Rmw = ^(l/A) 1 / 3 is the planet’s 
Hill radius. To avoid possible thermal shocks at the be¬ 
ginning of the simulation, we gradually inject the planet 
energy in the same way as we do for the planet mass (see 
Eq. EJ i.e., following Q+(£) = Q p sin 2 (7rt/10T p ), with T p 
the planet’s orbital period, and 

q+ _ { f(r)L p / (tt-Rcpd) if r = IV - f^| < R c pd 
1 0 otherwise, 

( 10 ) 

where f(r) is a Gaussian function used to smoothly in¬ 
ject the planet energy within the CPD. Its expression 
is f(r) = Aexp(—5 gt/R 0PD ), with A is a dimension¬ 
less normalization factor equal to 5. It has a FWHM of 
~ 0.75 Rcpd- 

It is worth noting that our model has inherent limita¬ 
tions due to its 2D geometry. For instance, the injected 
energy from the planet is only transported through the r- 
(j) plane of the CPD, neglecting the energy escape in the 
vertical direction. This could lead to an overestimation of 
the thermal energy deposited on the disk. 

2.3.3. Radiative cooling and disk spectrum 

The cooling term Q~ in Eq. m corresponds to the en¬ 
ergy flux radiated by the disk in the vertical direction. This 
quantity depends on whether the disk is optically thin or 
thick. Several processes could be responsible for the energy 
evacuation, e.g., for high temperatures ( rsj 10 4 K) Thomson 
scattering processes, free-free, bound-free transitions are 
dominant. For lower temperatures e.g., ~ 1 — 10 3 K, as 
in our case, Rosseland and Planck mean o pacities of dust 
and grain species sho uld be used (e.g., Bel l fc Lml (j 19941 ). 
ISemenov et all (|2003[ )). 

In our models, we calculate the optical depth (r = EE/2) 
assum ing the Rosseland mean opacity k from IBell fc Linl 
( 1994 ). In some regions of the disk, the gas density could 
be diluted (e.g., at the gap, and/or planet location), and 
become transparent or less opaque (r ^ 1). We use the 
iHubenvl (119901 ) prescription to calculate the effective opti¬ 
cal depth, 


7"eff 


3r 1 

~8 + — + 4t' 


( 11 ) 


The effective temperature (T e ff) and mid-plane temper¬ 
ature (T) are then related through 



( 12 ) 


The cooling rate per unit area due to radiation from the 
surface of the disk is calculated by integrating the total 
emitted radiative flux over all frequencies, 


poo 

Q~{r) = 2/ 4>^(T eff (r,0))dz/, (13) 

Jo 


where (r, 0)) is the local emergent flux of the disk 

surface. The factor 2 is needed because radiation escapes 
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from both sides (top and bottom) of the disk. For simplic¬ 
ity, we assume that the local emitted flux is given by the 
blackbody approximation, 

^(T) = 2^(T eff (r,0)), (14) 

where F>*,(T e ff(r, </>)) is the Planck function, and T e ff(r, 0) 
is the surface disk temperature. 

The spectrum of the disk is calculated by integrating the 
emitted flux 4v(T e ff(r, </>)) (Eq. [TTl) over the surface of the 
disk, 

L v = J <S>„(T eS (r,<j>))rdrd(t>. (15) 

Using Eq. (fl5j) we calculate the bolometric luminosity of 
the disk integrating over all the frequency domain, 



It corresponds to the electromagnetic energy per unit time 
radiated away in all wavelength (without including irradi¬ 
ation from the star). 

3. RESULTS 

In this section we present results of simulations for three 
different planet luminosities: L p = 10 -5 , 10 -4 , and 10 -3 
L 0 . We pay special attention to the gas properties (e.g., 
density, temperature) and to the spectral consequences of 
planet feedback for the highest planet luminosity, i.e. L p = 
1O- 3 L 0 . 

3.1. Density field 

We first examine how the density profile of the disk 
is affected by the planet’s radiative feedback. We com¬ 
pute the azimuthally-averaged density profile (E), given 

by (E) = ± / 0 27r E (r, <j))d<j>. In Figure U we compare two 
models: the first one has no feedback (L p = 0), while the 
second one assumes a planet luminosity of L p = 1O _3 L 0 . 
The density profile are shown at 300 planet orbits. At this 
stage, the gap carved by the planet is already formed and 
in a quasi steady state. 



Fig. 1.— Azimuthally-averaged density profiles without and with 
planet feedback (the planet’s luminosity is L p = 1O _3 L0). Both 
profiles are displayed after 300 orbits of the planet. Notice that the 
inner disk accumulates more material when feedback is activated. 

For this model, the effects of planet feedback are notice¬ 
able in a region of approximately 8 au radial extent about 


the planet. We notice that when the feedback is activated, 
the density profile increases in the inner disk (the disk re¬ 
gion inside the planet’s orbit), while the density decreases 
in the outer disk. Outside this region, there is no apparent 
difference if the planet emits energy or not. The observed 
changes in the density profile indicate that the planet’s lu¬ 
minosity enhances the ability of the disk to transport mass 
from the outer to the inner disk through the protoplane¬ 
tary gap. Although not shown here, we have checked that 
this effect is practically negligible for planet luminosities 
smaller than 1O -3 L 0 . 
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Fig. 2. — Density contours (cgs units) of the disk at 300 orbits. The 
top panel shows our results without planet feedback. The bottom 
panel is for a model with L p = 10 -3 L 0 . When feedback is activated, 
the inner disk is denser (bottom panel) than without feedback. There 
is an enhancement of the transport of matter from the outer to the 
inner disk when feedback is activated. From the figure one can also 
note that there is no matter accumulating at the circumplanetary 
region. 

In Figure [2j we compare density contours of the disk 
with and without feedback after 300 planet orbits. As 
previously shown by Figure [lj when feedback is activated 
the inner disk becomes slightly denser, suggesting that the 
flux of matter from the outer to the inner disk is enhanced 
by the planet feedback. 

It is worth mentioning that, when there is no feedback, 
the increment of the inner disk density is much smaller 
than in the case the feedback is activated (Fig. [2]). Recall 
that we use outflow (inner/outer) boundary conditions, 
hence the density increment is not due to material accu¬ 
mulating at the boundary, but rather a result of a change in 
the st elloc entric flux of matter stimulated near the planet 
(see § 13.51) . Notice also that in our models there is no 
material accumulating at the circumplanetary region. On 
the contrary, material is being slowly evacuated from the 
circumplanetary disk (as seen in Fig. [2j), but after 1000 or¬ 
bits there is still plenty of material producing an optically 
thick circumplanetary region. Correspondingly, the den¬ 
sity of the outer disk is slightly decreases when the planet 
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feedback is activated. We will show later that actually the 
disk’s stellocentric accretion rate is enhanced at the planet 
location, explaining this behaviour. 

3.2. Temperature field 



Fig. 3. — Azimuthally-averaged temperature profile without and 
with feedback (L p = 10 -3 L©). Both profiles are taken after 300 
orbits of the planet. 

In Figure [3] we show the azimuthally-averaged temper¬ 
ature profile (i.e., (T) = ^ T(r, of the disk for 

a model with a planet luminosity of 10 -3 L 0 and another 
without feedback, both after 300 orbits of the planet. We 
see that at the planet’s location the azimuthally-averaged 
temperature increases from about 50 K without feedback 
to nearly 80 K with feedback. Recall that these are 
azimuthally-averaged temperature profiles, and do not re¬ 
flect the large local variation shown below. 

In Figure 0] we compare the surface temperature of 
the disk after 300 orbits for the cases with L p = 0 and 
Lp = 10 -3 L 0 (note that a logarithmic scale is used, and 
that only the inner 20 au of the disk are shown). When 
no feedback is included, the maximum temperature oc¬ 
curs at the grid’s innermost radius. The maximum tem¬ 
perature reached at the planet’s location is about 150 K. 
When the feedback is activated, a hot spot forms within 
an au or so from the planet’s location, which is close to 
the size of the planet’s circumplanetary region (recall that 
in our feedback model, the energy released by the planet 
is injected in a region of area 7rR 0PD about the planet’s 
location, where Rcpd denotes the radius of the planet’s 
circumplanetary material, which is ~ 0.5 au). The maxi¬ 
mum temperature reached in the planet’s circumplanetary 
disk is about 1190 K. 

After 300 orbits of the planet, the aspect ratio at the 
planet location is about il p /r p ~ 0.2, therefore the pres¬ 
sure scale height of the disk at this position is H p ~ 2 au. 

3.3. Spectral signature 

We have shown in the previous subsection that the 
disk temperature near the planet’s location is strongly en¬ 
hanced by the inclusion of the planet feedback. This en¬ 
hancement (from ^ 100 K to ~ 1000 K) should cause a sig¬ 
nificant variation in the spectral emission of the disk in the 
vicinity of the planet. In Fig. Owe display after 300 orbits 
how the disk spectrum, calculated from Eq. 05J changes 
for different planet luminosities (L p = 10 , L p = 10 -4 , 

Lp = 10 -3 , and Lp = 0 L 0 ). We see that the spectrum 




Fig. 4. — Effective temperature of the disk after 300 orbits of the 
planet. The top panel shows our model without planet feedback, 
the bottom panel is for our model with L p = 1O _3 L 0 . Without 
feedback, the circumplanetary region reaches temperatures of 160 K, 
while when the planet luminosity is included this area reaches a peak 
temperature of about 1190 K. 

is dramatically modified when the feedback is higher than 
L p > 1O- 4 L 0 . 


100 orbits 



Fig. 5.— Spectrum of the disk after 300 orbits without and with 
planet feedback, for L p = 10 -5 , L p = 10 -4 and L p = 10 -3 solar 
luminosities. 

At the planet’s location the disk temperature without 
feedback is about 150K. When the feedback is included, 
the temperature derived from the energy equation in the 
vicinity of the planet peaks at 1100 K0 for planet luminosi¬ 
ties L p > 1O -3 L 0 . This increase in temperature produces 

2 We point out that the high te mpera tures we observe in our nu¬ 
merical models are also recovered bv IZhul (120151 k In his 1-dimensional 
accretion disk models the temperature down at the atmosphere of the 
planet (r ~ Rj) reaches ~ 2000 K, without invoking any energy input 
from the planet itself. Notice that those scales cannot be resolved 
with our numerical scheme, so a direct comparison is not possible. 
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a corresponding increase in the peak flux density by a fac¬ 
tor of ~ 100, as seen in Fig. [5] Notice also that without 
feedback the spectrum of the disk peaks at 19.3/im, and 
that for a planet luminosity of L p = 10 -3 L 0 , the spectrum 
peaks at 2.5/im. 

Using Eq. [16j we can integrate the spectrum from Fig. [5] 
to obtain the bolometric luminosity of the entire dislfl In 
the simulation without feedback, the bolometric luminos¬ 
ity is Lb = 5.54 x 1O _2 L 0 , while in the case when we 
assume a planet luminosity of L p = 10 -3 L 0 , the bolomet¬ 
ric luminosity increases to Lb = O.98L 0 . In other words, 
the radiative feedback from the planet results in a disk 
luminosity increased by a factor of 17.8. 

To better appreciate the differences between the cases 
L p = 0 and L p = 1O -3 L 0 , we plot in Fig. [6] (log scale) 
the ratio between the bolometric emission per unit area 
(given by f^° <F(T(r, </>))dA = crT 4 ) for L p = 10 -3 L 0 and 
the emission when L p = 0, after 300 orbits of the planet. 
From Fig.[6]we notice that when the feedback is activated, 
the blackbody emissions (crT 4 ) at the vicinity of the planet, 
inside a radius of about ~ 5 au, increases by up to 3.6 or¬ 
ders of magnitude over a situation without feedback. Also, 
we notice that along the orbital path of the planet around 
the central star, the feedback leaves a ‘track’ of enhanced 
disk surface brightness by about one order of magnitude. 
This track indicates that the radiative feedback from the 
luminosity of the planet induces net heating of the gas not 
only in the close vicinity of the planet, but also along the 
planet’s trajectory. 



x[AU] 


Fig. 6.— Map ratio (log 10 scale) of the bolometric emission per 
unit area between our model with L p = 10 -3 L 0 and our model with 
L p = 0, at 300 orbits. When the feedback is activated, the bolometric 
surface brightness is about 10 3-6 times higher in the circumplanetary 
region, and there is also an excess of emission along the planet’s orbit. 

In Figure [3 we display at 300 planet orbits the black- 
body radiation B\(T(r,4>)) for A = 2.5/im for the cases 
without and with feedback (L p = 1CR 3 L 0 ). The feedback 
dramatically increases the blackbody radiation in the cir¬ 
cumplanetary region. We point out that there could ap¬ 
pear to be an energy conservation problem here, because 
when there is no feedback we get a luminosity output from 
the entire disk of about Lb = 5.5 x io -2 L©, while adding a 
relatively small source of L p = 1CR 3 L 0 , we obtain an out- 

3 Notice that this bolometric luminosity does not include the emis¬ 
sion from the star. Moreover, the quoted values depend on the radial 
extent of the disk, most importantly on its inner radius, which we 
set for numerical convenience to 1 au. 


put luminosity of O.98L 0 . In the next subsection, we run 
several tests in order to check this result, concluding that 
there is no inconsistency with energy conservation. The 
additional energy comes from an increase on the stellocen- 
tric accretion rate through the di sk w hen the feedback is 
activated, as will be explained in £13.51 below. 



x[AU] 


Fig. 7. — Local blackbody flux per unit area B\(T{r, < f >)) °f the 
disk according to Eq. [Til taking A = 2.5 fim. The map is in log scale 
and the units are erg s -1 cm -3 . The upper panel corresponds to a 
model without feedback, the bottom panel is for L p = 10 _3 Lq. The 
highest surface brightness comes from the planet’s circumplanetary 
region. 


3.4. Feedback behaviour for other viscosity prescriptions 

In this section we examine to what extent the impact of 
the planet’s radiative feedback depends on the assumption 
on the viscosity prescription. For this purpose, we have 
carried out a number of simulations with different viscosity 
prescriptions. 

We first employed alpha viscous disk models. The alpha 
disk model, as it is implemented in the public release of 
FARGO-AD, assumes a kinematic viscosity v = a(c%/Q) 
where the brackets stand for the azimuthal average. Dif¬ 
ferently said, the kinematic viscosity is purely radial. We 
have carried out simulations with « = 4x ICR 4 , 4 x ICR 3 
and 4 x ICR 2 , and found nearly identical peak tempera¬ 
ture in the planet’s circumplanetary region. Our results 
are listed in Table [lj along with disk’s accretion rates at 
the planet location, which we will detail in the next sec¬ 
tion. 

Lastly, we ran an inviscid model, i.e., a model without 
shear viscosity (v = 0; note, however, that this model still 
includes artificial viscous heating via a bulk viscosity, as 
in all the m odels p resented in this paper - see model de¬ 
scribed in § 12.3.1) ) . In that case, we find again a peak 
temperature near the planet of 1166 K, which is very sim¬ 
ilar to the cases presented before with different viscosity 
prescriptions. All these numerical experiments show that 
the peak temperature that we find with planet feedback 
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activated is basically independent of the disk’s viscosity. 


Viscosity model 

Loeafc [K] 

T^ak [K] 

Mon[M 0 /yr] 

Moff[M 0 /yr] 

a = 4 x nr 4 

1170 

125 

3.7 x 10“ 6 

2 x in 8 

a = 4x nr 3 

1190 

200 

2.6 x 10“ 6 

1.2 x 10“ 7 

a = 4x 10 -2 

1200 

276 

3.8 x 10“ 6 

2.3 x 10“ 7 

v = 0 

1166 

100 

1 x 10“ 6 

1 x 10“ 8 


TABLE 1 

Peak temperature and disk accretion rate values at the 

CIRCUMPLANETARY REGION FOR DIFFERENT VISCOSITY 
PRESCRIPTIONS, WITH FEEDBACK (COLUMNS 1 AND 3 ON THE 
RIGHT-HAND SIDE, L p = 1O _3 L 0 ) AND WITHOUT FEEDBACK 
(COLUMNS 2 AND 4). 


3.5. Accretion rate 

We have shown in S l3.3l that the bolometric luminosity of 
the disk increases by a factor of ~ 20 with planet feedback 
and a planet luminosity of L p = 1O _3 L 0 , compared to a 
disk without planet feedback. This excess of energy mostly 
arises from the planet’s circumplanetary region and from a 
narrow region about the planet’s orbital radius (see Figure 
Ej). The luminosity can be related to an accretion rate, an 
excess of luminosity should therefore have a corresponding 
increase in the disk’s accretion rate. 



Fig. 8. — Disk accretion rate at the planet’s orbital radius, defined 
by Eq. (Q2J, without and with (L p = 10 3 L 0 ) planet feedback. 
With feedback the accretion rate at the planet’s orbital radius is 
enhanced by an order of magnitude. The peak of the accretion rate 
is reached near ~ 40 planet orbits, at this stage the temperature 
reaches a maximum near the planet vicinity. 

We calculate the azimuthally-integrated accretion rate 
of the disk at the planet location as 

M = J v r {r v ,<t>)T,{r p ,(j))r v d<j), (17) 

where v r is the radial velocity of the gas. In Figure [8] we 
display the time evolution of M for two models: one with 
feedback (L p = 10 -3 ) and one without feedback (L p = 0), 
assuming an alpha disk with a = 4 x 10 -3 (our fiducial 
model). We see that, when the feedback is activated, the 
disk’s accretion rate is enhanced by an order of magnitude. 
The accretion rate reaches similar values for different vis¬ 
cosity models, see Table HJ 

We notice from Figure |T] that the density average at 
the planet location did not differ much between the cases 
with and without feedback, hence, in order to have an en¬ 
hancement of the accretion rate at that region, it must 
be the radial velocity that is enhanced. This is exactly 
what we see in Figure [9j where the radial velocity is en¬ 
hanced by a factor ~ 4 when the feedback is activated. It 


is not expected to have a similar enhancement factor as 
the observed for the accretion rate mentioned above (i.e., 
one order of magnitude) because Figures [Hand [9] show an 
azimuthally-average radial velocity, rather than local val¬ 
ues at the planet location. To compute M(r p ) as indicated 
in the precedent paragraph we use local values (where, for 
instance, the density takes larger values at r p than the 
azimuthally averaged value), then we integrate over the 
azimuthal coordinate. 



Fig. 9. — Azimuthal average of the radial velocity after 600 orbits. 
When the feedback is activated, the radial velocity is enhanced by a 
factor ~ 4 at the planet location (10 AU). 

At the end of we discuss the fact that our simula¬ 

tions show, for instance, that when we use a local energy 
input of L p = 1O _3 L 0 , we obtain an energy output of 
Lb ^ L 0 . It could be misinterpreted as a non-conservative 
evolution of the energy. To show that it is not the case, 
we argue that the origin of the extra energy comes from 
an increase in the disk’s accretion rate stimulated by the 
feedback of the planet, and that this enhancement is inde¬ 
pendent of the viscosity prescription. 

When radiation feedback is active, the luminosity of the 
disk at the location of the planet is proportional to the 
stellocentric disk accretion rate (which attains a maximum 
value), i.e., L acc oc M at R p . For instance, taking the 
model with a = 4 x 10 —3 , the accretion rate at 300 or¬ 
bits is M on ~ 10 -6 M o yr -1 with feedback, and M off = 
6 x lO -8 M 0 yr -1 without (see Figure [8]). From this, we 
find that the accretion luminosity of the disk should be in¬ 
creased by a factor of roughly L°£ c /L°^ c = M on /M off ~ 17 
by the inclusion of planet feedback. This is indeed in good 
agreement with the increase in bolometric luminosity cal¬ 
culated in § 13.31 where the luminosity increase is found to 
be « 17.8. Our results clearly show that there is no im¬ 
balance in the energy bill during the simulation, and that 
the extra energy is the direct product of an enhancement 
in the disk accretion rate, which in turn was stimulated by 
the planet feedback. 

Now we briefly describe the effect on the disk dynam¬ 
ics when the radiative feedback of the planet is included. 
From a dynamical point of view, the radial component 
of the Navier-Stokes equation includes the radial pressure 
gradient of the fluid (the term (1 /YAjdP/dr). A local en¬ 
hancement of this term at the planet location should lead 
to a local increase in the radial velocity (understood in ab¬ 
solute value). In Figure [101 we compare the azimuthally- 
averaged radial acceleration due to the gradient pressure 
when the feedback is turned on (L p = 1O _3 L 0 ) and off, 
after 100 planet orbits. We show for comparison the radial 
gravitational acceleration due to the central star. We see 
that the pressure gradient is particularly strong near the 
planet’s location, as expected. Moreover, with feedback 






























Fig. 10.— Azimuthal average of the radial acceleration due to the 
pressure work done by the disk without and with feedback. The 
gravitational radial acceleration is displayed for comparison. 

the radial acceleration is more negative, which acts to in¬ 
crease the accretion rate through the disk near the planet’s 
location, as seen before. 

Similar conclusions were obtained bv IQwen ( 2014 ). who 
construct a ID radial ’transition’ disk model including 
feedback from the protoplanet acting over the gas and dust 
components of the disk. They include the radial and az¬ 
imuthal accelerations of the dust particles produced by 
radiative feedback from the planet in a secular model of a 
proto planetary disk with an embedded accreting planet 
(e.g., IClarke fc Pringlel (|1988h ). finding that the accre¬ 
tion rates observed in transition disks are better explained 
wh en t he ra diative feedback of the planet is included. The 
IQwenl (120141 ) model assumes azimuthal symmetry, in which 
the radiation feedback and the disk-planet interaction are 
treated as average quantities. By contrast, our model is 
non-axisymmetric (2D), and we include a non-stationary 
energy equation to introduce the radiative feedback. 

It is important to clarify that this accretion enhancement 
should lead to a transient situation explained as this: as 
shown above, the local feedback promotes the flux of mat¬ 
ter in a region close to the circumplanetary disk. As seen 
in Figure [2j gas density slowly diminishes at the CPD. If 
this continues long enough, the CPD will be depleted be¬ 
coming optically thin. If this happens, radiation from the 
planet will no longer interact with the gas, escaping im¬ 
mediately from the disk without heating it. Therefore, the 
radiative feedback will no longer affect the CPD, muting 
its effect on the gas dynamics. Consequently, accretion 
enhancement will be no more, as well as the produced ex¬ 
tra luminosity. Once the CPD start again to accumulate 
material and became anew optically thick, radiation feed¬ 
back will interact one more time with the gas, modifying 
its dynamics as before. From our calculations, to deplete 
the CPD with a planet feedback of 10- 3 L 0 located 
at 10 au, takes more than 10 thousand years, but further 
investigations are needed to clarify this point. 

To summarize, our findings show that a relatively small 
source of heating coming from the planet’s circumplane¬ 
tary region will induce enhanced disk accretion rates near 
or around the planet’s orbital radius. This enhancement 
facilitates the extraction of gravitational energy, which lo¬ 
cally heats the gas, enhancing the luminosity in the cir¬ 
cumplanetary region. This local increase in the disk tem¬ 
perature then leads to an increase in the accretion rate, 
increasing even more the temperature (and thus the lu¬ 
minosity). This effect results in a positive feedback until 
a thermal balance is reached where the dissipation heat 
rate equals the radiative cooling rate (Q + = Q~ with our 
notations). 


These accretion luminosities should be detectable inside 
gaps. It should be noticed that it is not the intrinsic lumi¬ 
nosity of the planet, but a local gas luminosity stimulated 
by the action of the planet feedback. In the next sub¬ 
section, we present results of hydrodynamical simulations 
dedicated to interpret the observations of a protoplanet 
candidate around star HD 100546. 

3.6. HD 100546 simulation 

A protoplanet injecting extra local heating in the disk via 
the feedback mechanism explained in the previous sections 
could explain the bright compact emission d etecte d in L' 
band with NACO in the disk of HD 100546 (IQuanz et al.l 
12013a ). HD 100546 is a Herbig Ae/Be star which harb ours 
a prot o planet candi d ate orb iting at 68 au ([Currie et al.l 
([2014D , IQuanz et al.l (|2013b[) ). Motivated by these obser¬ 
vations, we have carried out a set of simulations for HD 
100546. We compare the results of various accretion lumi¬ 
nosities with the luminosity of the observed planet candi¬ 
date. 

Based on interferometric d ata u sin g AMBER /VLTI and 
photometric observations, ITatulli et al.l ( 201 1) proposed 
a disk model for the circumstellar environment of HD 
100546. Using this model we adopt a disk scale height pro¬ 
file H{r ) = 12 (r/100 au) 1,1 au, and a surface density pro¬ 
file £(r) oc r -1 as initial conditions for our simulations. We 
tailor the disk mass in order to guarantee that the circum¬ 
planetary region remains optically thick. We thus set the 
disk gas mass to Md = 15 x 1O _2 M 0 (while ITatulli et al.l 
([2011 ) reported 5 x 1O _2 M 0 for gas). From [Panic et al.l 
(2010), we assume a central star of mass M * = 2.5M 0 
(which defines our unit of mass), implying a planet’s or¬ 
bital period of 354.6 yr. Neither irradiation nor flux dif¬ 
fusion from the cent ral star are taken into account in our 
simulation. In §EU we justify this approximation for the 
irradiation through a simple analytical calculation. As in 
our fiducial model, the disk is treated as an alpha viscous 
disk with a = 4 x 10 -3 . 

The planet is located at r p = 68 au, which is taken as the 
code’s unit of length. For numerical convenience, the grid 
now extends from 10 to 200 au along the radial direction, 
which corresponds to radii ranging from 0.15 to 3.0 in code 
units. 

In forme d by our dust model for HD 100546 (see 
Sec. 13.6.11 below), we changed the opacity prescription in 
our feedback calculation taking into account an average 
opacity for a mix of dust species given by k = 130 cm 2 /g. 
This is more adequate for transition disks such as HD 
100546. 

We carried out a set of simulations with various planet 
luminosities: L p = {0; 2.5; 5} x 1O _4 L 0 . The observed 
emission in L' band is thought to correspond to a p roto¬ 
planet with a mass in between 1 and 8 Mj ([Quanz et al.l 
I2013al) . We fix the planet mass to M p = 5Mj. Our re¬ 
sults of simulations are shown after 600 planet orbits, or 
24.2 x 10 3 yr (as in the previous sections, the planet’s mass 
reaches its imposed value over 10 orbital periods). 

We show in Figure fTTI contours of the disk temperature 
for the different planet luminosities mentioned above. The 
coordinates in the figures are in code units, therefore the 
planet (at 68 au) appears at r = 1. 

The top panel corresponds to a simulation without feed¬ 
back. In this case, the hottest temperatures come from 
the innermost region of the disk (200 K). The circum¬ 
planetary region reaches a maximum temperature of about 
80 K. The middle panel shows a model with feedback and 
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a planet luminosity L p = 2.5 x 1O _4 L 0 . In this case, the 
proto-planet region peaks at 270K. In the lower panel, 
L p = 5 x 1O“ 4 L 0 , which gives a peak temperature (at the 
planet location) of 325K. 


600 orbits - Temperature [K] ( L p =0) 



-3-2-10 1 2 3 

x [code units] 

600 orbits - Temperature [K] ( L p =2.5 xlO“ 4 L 0 ) 



-3-2-10 1 2 3 

x [code units] 

600 orbits - Temperature [K] (L p =5 xl(T 4 L 0 ) 



-3-2-10 1 2 3 

x [code units] 
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Fig. 11. — Contours of the disk temperature after 600 planet orbits. 
From top to bottom: L p = 0, L p = 2.5 X 10 _4 Lq, and L p = 
5 X 10- 4 Lq. The peak temperature near the planet is 80 K, 270 K 
and 325 K, respectively. 


By integrating the spectrum over the whole disk surface 
and over the entire wavelength range (Eq. fl6|) . we obtain 
the bolometric luminosity of the disk. Without feedback, 
it yields = 4.2L 0 . For L p = 2.5 x 1O _4 L 0 , we ob¬ 
tain a disk luminosity of about L b = 4.7 L 0 , resulting 
in a disk 1.1 times brighter than without feedback (i.e., 
Ton/Toff = 1.1). This can be explained by the fact that 
the addition of a small source of heating near the planet 
(associated with L p = 2.5 x 10- 4 L g ) locally enhances the 
disk accretion rate at the planet location, and results in 
a positive feedback able to produce a disk that is ^ 10% 
brighter than without feedback. We point out that the 
disk luminosity without feedback is 4.3 solar luminosities 
(for a disk extending from 10 to 200 au), which is consis¬ 
tent with the bolometric luminosity of ^10 T 0 reported 
by Benisty et al. (2010) (for a disk from 10 to 500 au). 


Recall that we neglect irradiation from the star, and that 
modelling a smaller than observed disk is likely to explain 
the factor ~ 2 discrepancy. 

3.6.1. Radiative transfer 

In order to compar e our HD 100546 simu lation with the 
L' observations from IQuanz et al.l ( 2013b ). we input the 
hydrodynamical gas density and gas temperature fields 
into t he RADMC3DB radiative transfer code (version 
0.38, (jDullemond et al.l 120141 )). At infrared wavelengths 
the emission is dominated by thermal and scattered emis¬ 
sion from micron-sized dust. Small grains should be well 
mixed, hence we assume that the dust distribution follows 
the same density field as the gas in the simulation. These 
grains should account for the bulk of the L' emission. 

Our model dust distribution consists of a mix of 2 com¬ 
mon s pecies: amorphous carbons and astronomical sili¬ 
cates (Dr aine fe Le e 1984). As i nformed by p r eviou s SED 
modeling (jBenistv et al.l ( 201(1 ). ITatulli et all (|2011[ )). our 
grain size distribution follows a power-law with exponent 
—3.5, and particle sizes ranging from 0.05 to 1000 /im. We 
used Mie theory (homogeneous spheres) to compute the 
dust opacities for anisot ropic scattering. The o ptical con¬ 
stants were taken from ILi fe Gre enberg (119971) for am or- 
phous carbon grains, and from IlDraine fe Lea (|1984f ) for 
silicates. The intrinsic densities of the grains are 2 g cm -3 
and 4 g cm -3 for amorphous carbons and silicates, respec¬ 
tively. We also assume Tgas — Tdust- 

The 2D surface densities for each dust species are ex¬ 
tended vertically following hydrostatic equilibrium. This 
produces a 3D volume with a Gaussian profile in the ver¬ 
tical direction. The disk scale-height is assumed to be the 
same for both dust species and it is set by the temperature 

field through the relation H (r) = yJ^RT{r) r /^Kep, where 

VKep is the Keplerian velocity. For simplicity the tempera¬ 
ture was assumed constant in the vertical direction. A bet¬ 
ter vertical description of the temperature would require a 
full account of stellar irradiation and stellar flux diffusion 
in the simulations. The final image is produced performing 
a second order volume ray-tracing method. The results 
of our radiative transfer calculation for L' are shown in 

Fig.m 

IQuanz et all ( 2013b ) calculated the luminosity of the hot 
spot around the companion candidate HD 100546b, assum¬ 
ing a point source embedded in a region 0.1 arcsec in size 
(^ 10 au) located at 68 au from the central star, reporting 
an apparent L' magnitude of 13.2, which translates into 
a flux density of 1.3 mJy. In order to compare with our 
planet feedback predictions, we computed the emergent 
flux from the vicinity of the protoplanet in the same way. 

The simulation including an embedded protoplanet with 
an accretion luminosity of T p = 2.5 x 1O -4 L 0 yields an 
emerging flux from the circumplanetary vicinity of 0.7 mJy 
in L '. On the other hand, an accretion luminosity of 
L p = 5 x 10 -4 T 0 already produces a flux density of ^3 Jy, 
abo ut two times hig her than the flux levels reported by 
IQuanz et al.l (|2013bh . 

A key result from our model is that the emitted luminos¬ 
ity from the protoplanet region depends on the prescribed 
feedback and not on the viscosity prescription. Therefore, 
from our simulations for HD 100546b, we conclude that 
the proto-planet candidate is compatible with an accret- 

4 http://www. ita.uni-heidelberg.de/~dullemond/software/radmc- 
3d/ 
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Fig. 12.— Radiative transfer predictions in L' band (3.8 jum), for a disk inclination angle of 47 deg, calculated on the density and temperature 
fields shown in Fig. 1111 The left panel shows the results with radiative feedback (L p = 2.5 x 10 -4 L 0 ), the right panel shows the case without 
feedbac k. The radiative feedback results in a circumplanetary hotspot reaching ~1 mJy in L ', close to the 1.3 mjy observed bv lQuanz et al.l 
(l2013bT ). The inclusion of the inner disk (with radius <10 au) would result in a bright stellar point source, which could efficiently be cancelled 
by Angular Differential Imaging (ADI). The images have been convolved with a 0.1 arcsec Gaussian PSF. Colour stretch is logarithmic. Images 
are in units of ergcm _2 s —1 str -1 . 


ing 5 mass Jupiter-like object, with a planet luminosity of 
about ~ 2.5 x 1O _4 L 0 . 

4. DISCUSSION 

4.1. Irradiation from the central star 

Here we briefly discuss the possible effect of the irradia¬ 
tion from the central star on the disk. The flux penetrating 
the_surface of the disk can be calculated as (see for instance 
IFrnhlichl (12001 ') 


F irr = (1 -/3)f^ cosip, (18) 

where /3 is the albedo (reflection coefficient), and (p is the 
angle formed by the incident radiation and the normal to 
the surface. It can be shown that cos(^ cx dH/dr — H/r, 
and using F irr = <7l)< we can compute the irradiation 
temperature T- lvv at the surface of the disk as 




I/* fH\ fd\aH 
4nr 2 cr \r J \ d\nr 


l) (1 - P) • 


(19) 


From analytic (Chi ang fc Golcxeichl 119971) and numer¬ 
i cal s olutions (D’Ale ssio et al.l ( 19981 ): ibullemond et al.l 
(J2002D ). it can be shown that the height H of the disk 
when stellar radiation is taken into account has a power- 
law dependence with radius, H ex r^, with / ~ 1.3 — 1.5. 
Using these values, we obtain that (d In H/d In r — 1) varies 
slightly between 0.3 to 0.5. The aspect ratio ( H/r ) can 
reach values from 10 -4 to 1CT 1 (e.g., IBell et all ([19971 )). 
Even taking values that give the maximum temperature 
from the above equation i.e., dlnH/dlnr — 1 = 0.5; 
H/r = 10 -1 ; and /3 = 0 (all rad iatio n penetr ates), and as¬ 
suming a central star of 50 L 0 ([Panic et al.l ( 2010 ) quotes 
26 L 0 for HD 100546), the disk temperature from irradia¬ 
tion at r = 68 au reaches ^ 60 K, while from viscous dis¬ 
sipation plus radiative feedback the temperature reaches 
values > 100 K depending on the planet luminosity L p . 

For the simulations with a solar-type star presented in 
§E3H3jg we obtain a surface temperature of about ^ 55 K 
at r = 10au (planet location). Therefore, the heating of 


the disk by the central star is, in the most extreme case, 
of the same order as the temperature obtained by viscous 
dissipation when there is no feedback (T ~ 55K), and one 
order of magnitude bellow when the feedback is activated 
T ~ 1000K (see Figured]). 

We conclude from this that the surface temperature in 
our 2D disk model is not expected to be affected by the 
irradiation energy of the central star, and that this effect 
can be neglected from our calculations. 

4.2. Application to massive black hole binaries 

Even though in this paper we focus on proto-planetary 
discs, a very similar physical set-up is that of unequal-mass 
massive black hole binaries. Such binaries are expected to 
form after galaxy mergers, once the central massive black 
hole of each galaxy migrates to the ce nter of the new 
system due to dynamical friction (e.g., Begelman et al.l 
([19800 ). During the merger, large quantities of gas get 
funnelled to the inner region of the new galaxy, where they 
form_a nuclea r gas d isk that interacts with the binary (e.g., 
IMaver et al.l ([20071) ). If the masses of the black holes are 
dissimilar, the secondary will orbit the primary, while both 
are embedded in the gaseous disk, in a situation very much 
resembling the star - planet-disk systems we model in this 
paper fe.g. JArmitage fc Nataraianl ( 2002 )). especially con¬ 
sidering that both planets and embedded black holes are 
expected to produce luminosity by accretion. 

Binaries with separations of parsecs or smaller are not 
directly resolvable, so we need to rely on indirect meth¬ 
ods to identify them. One such method is based on the 
assumption that the gaseous disk around the primary ra¬ 
diates as a multi-color blackbody. If the secondary pro¬ 
duces a gap in the disk, then the spectrum will show a dip 
at the wavelength range associated w it h the temperature 
of the missing gas ([Chang et al.l ([20101 ): iGultekin fc Millerl 
([2012D ). Extrapolating our results to that regime, it is 
clear that the disk spectra will also be modified by the 
heating up of the gas surrounding the secondary black hole, 
likely producing a larger dip and a concurrent increase of 
the shorter wavelength emission. A more detailed analysis 
of this situation is deferred to a follow-up paper. 
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5. SUMMARY 

In this paper we present 2D hydro dynamical simulations 
of the interaction between a Jupiter-mass planet and its 
parent protoplanetary disk, including for the first time the 
effect of the planet’s radiative feedback onto the disk (en¬ 
ergy released by the planet as it assembles its mass). In 
this first work, we assume that the energy released by the 
planet to the disk is constant over the duration of our 
simulations, a few hundred planet orbits typically. We 
carried out various simulations taking realistic parameters 
for protoplanetary disks that could be identified wit h actu ¬ 
ally ob s erved system s , like e .g., HD 100546 (|Currie et al.l 
(2014), IQuanz et al.l (|2013afl h varying the luminosity of 
the planet. 

We find that planet luminosities below L p < 10 -4 L© 
barely modify the disk response to the planet: the gap 
formed by the planet does not show any noticeable modifi¬ 
cation compared to without feedback and the emitted spec¬ 
trum of the disk remains practically unchanged. However, 
planets with L p > 10 _4 L© introduce significant changes: 
the additional energy input from the planet heats up the 
disk, locally enhancing the local disk accretion rate (under¬ 
stood into the primary), increasing even more the temper¬ 
ature at this location. This mechanism results in a positive 
feedback magnifying the energy output until a thermal bal¬ 
ance is reached. 

Figure [8] shows the increase in the disk accretion rate 
in the planet vicinity. The spectrum and luminosity of 
the disk are modified, shifting the emissions to higher am¬ 
plitudes and shorter wavelengths. For instance, assuming 
10 -3 L© for the planet feedback, the accretion rate at the 
planet vicinity increases from 6 x 10 -8 (without feedback) 
to 10 -6 M©yr _1 , and translates into an increase in the 
disk luminosity from Ldisk = 5.5 x lO -2 L© (without feed¬ 
back, Apeak — 19.3/1777/) to T^isk — 0.9 L© (A pe ak — 2.5/7777/). 


We find that our results do not depend on the viscos¬ 
ity prescription. Thus, in our models, setting the planet 
luminosity will fix the accretion rate, temperature and lu¬ 
minosity in a region close to the planet. Our results imply 
that observations of protoplanetary disks where planets are 
form could reveal the accretion process onto them, with¬ 
out much interference from nuance parameters such as the 
viscosity. 

Finally, we build a mode l for the system around HD 
100546 (Qua nz~et al.lI2013alh reproducing quite well the 
observed flux density in L' band of 1.3 mJy from the region 
around the accreting candidate planet HD 100546b. Our 
model indicates that this system contains a forming planet 
with a luminosity of rsj 10- 4 l©. 
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